if (interactive()) {
  source("./pipe.R")
  input = get_input(get_current_path())
  output = get_output(get_current_path())
}

paged_table(data.frame(value= unlist(params)))

1 Data

library(TCGAbiolinks)

## Warning: multiple methods tables found for 'sort'
## Warning: multiple methods tables found for 'sort'

library(biomaRt)
library(kableExtra)
library(rmarkdown)
library(readxl)
source(input$functions)
source(input$fc_params)
suppressWarnings(library(GSEABase))

## Loading required package: BiocGenerics

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union

## The following objects are masked from 'package:igraph':
## 
##     normalize, path, union

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter,
##     Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff, table, tapply, union, unique, unsplit, which.max, which.min

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and
##     for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:rlang':
## 
##     exprs

## Loading required package: annotate

## Loading required package: AnnotationDbi

## Loading required package: stats4

## Loading required package: IRanges

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:clusterProfiler':
## 
##     rename

## The following object is masked from 'package:plyr':
## 
##     rename

## The following objects are masked from 'package:dplyr':
## 
##     first, rename

## The following object is masked from 'package:tidyr':
## 
##     expand

## The following objects are masked from 'package:Matrix':
## 
##     expand, unname

## The following object is masked from 'package:utils':
## 
##     findMatches

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:clusterProfiler':
## 
##     slice

## The following object is masked from 'package:plyr':
## 
##     desc

## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice

## The following object is masked from 'package:purrr':
## 
##     reduce

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:clusterProfiler':
## 
##     select

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: XML

## Loading required package: graph

## 
## Attaching package: 'graph'

## The following object is masked from 'package:XML':
## 
##     addNode

## The following object is masked from 'package:plyr':
## 
##     join

## The following object is masked from 'package:stringr':
## 
##     boundary

## The following objects are masked from 'package:igraph':
## 
##     degree, edges, intersection

library(grid)
library(ComplexHeatmap)

## ========================================
## ComplexHeatmap version 2.10.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##   genomic data. Bioinformatics 2016.
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================

library(ggpubr)
library(stringr)
library(data.table)
library(ggstatsplot)

## You can cite this package as:
##      Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
##      Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167

# load hpv signature genes
avg_log2FC_cutoff = log2(signature_fc)
fdr_cutoff = signature_fdr


opscc_deg = read.table(file =  input$opscc_deg,row.names = 1,sep = ",",header = T)
hmsc_deg = read.table(file = input$hmsc_deg,row.names = 1,sep = ",",header = T)

hmsc_signatue = rownames(hmsc_deg[hmsc_deg$fdr<fdr_cutoff & hmsc_deg$avg_log2FC>avg_log2FC_cutoff,])


opscc_signature = opscc_deg %>% filter(avg_log2FC > avg_log2FC_cutoff &
                                  fdr < fdr_cutoff) %>% rownames()

message("opscc_signature:")

## opscc_signature:

print(opscc_signature)

## [1] "HPV16-E1" "HPV16-E7" "HPV16-E5" "FBL"      "DST"      "ALDH3A1"  "IGFBP2"   "KRT15"    "SERTAD1"

message("hmsc_signatue:")

## hmsc_signatue:

print(hmsc_signatue)

##   [1] "ANLN"          "TUBB6"         "AXL"           "IFT122"        "GFM1"          "SEPT2"         "PBRM1"         "INSIG1"       
##   [9] "CAPG"          "POP7"          "BHLHE40"       "MCM5"          "TRAPPC3"       "MRPS6"         "GINS2"         "NARS"         
##  [17] "LSM5"          "GUSB"          "EZR"           "FZD6"          "STAT2"         "TMPO"          "SLC35F5"       "MVB12A"       
##  [25] "CDC42EP3"      "COG4"          "MRPL20"        "GSDMD"         "VCL"           "MT1E"          "TUBA1C"        "TXN"          
##  [33] "ATXN10"        "USP39"         "MSH6"          "SDC2"          "OAT"           "SCAP"          "PRELID3B"      "KNTC1"        
##  [41] "ARF6"          "SMC5"          "DCTN2"         "AAMP"          "TBCB"          "KLF10"         "TPP1"          "IFT52"        
##  [49] "ETS2"          "EWSR1"         "NOP10"         "KRCC1"         "ASL"           "IGSF9"         "PHF5A"         "XRCC5"        
##  [57] "UBE2T"         "IK"            "AC007250.3"    "YES1"          "SNRPD1"        "C9orf3"        "FDFT1"         "TPM3"         
##  [65] "LASP1"         "SNRPC"         "PACSIN2"       "MTHFD2"        "PMVK"          "OBP2B"         "AKR1B1"        "SYPL1"        
##  [73] "LPAR2"         "PHB"           "THOC6"         "TMX1"          "SDHA"          "NUBP2"         "NDUFS7"        "TAOK3"        
##  [81] "MRPS9"         "PTRF"          "PSME2"         "MED15"         "PES1"          "RP11-19E11.1"  "RP5-1102E8.2"  "SEH1L"        
##  [89] "MCM3"          "POLR2L"        "DNAJC8"        "COPS6"         "PRKCDBP"       "KRT16"         "RNASEH2A"      "PDIA6"        
##  [97] "RAB11FIP1"     "RFC3"          "POLD2"         "AC079337.1"    "NF2"           "BABAM1"        "TPGS2"         "LOX"          
## [105] "SERINC3"       "NDUFB8"        "CLIC4"         "SELT"          "GRB2"          "HELLS"         "CXCR4"         "PRDX5"        
## [113] "FLII"          "EDIL3"         "SF3A3"         "PSMB1"         "FARSA"         "PREB"          "POMP"          "MPPED2"       
## [121] "CKB"           "MCOLN1"        "MCM4"          "PCNT"          "BTF3L4"        "DDX1"          "BRE"           "GPN1"         
## [129] "H2AFZ"         "KIF22"         "PSMB4"         "TMEM165"       "RP11-45P22.1"  "ECH1"          "NUCB2"         "NUDCD3"       
## [137] "C19orf10"      "VPS4B"         "MGME1"         "TCF19"         "KRT5"          "ATP6V0B"       "GTF2I"         "TSC2"         
## [145] "CNOT7"         "UFD1L"         "SMARCA5"       "TCERG1"        "PLAUR"         "SSRP1"         "SMC3"          "RASSF1"       
## [153] "TYMS"          "GPC4"          "CNPPD1"        "PLXDC2"        "UBXN11"        "RAD23A"        "USP4"          "C16orf13"     
## [161] "CTNNAL1"       "LRRCC1"        "ADNP"          "SDHB"          "UQCRQ"         "MTHFD1"        "MRPL21"        "RANGAP1"      
## [169] "PSMA5"         "BRIX1"         "RP11-665C16.9" "FDPS"          "FHL2"          "MYC"           "KLHDC3"        "DNAJB6"       
## [177] "RP5-1085F17.3" "DCTPP1"        "RBM8A"         "UBE2L3"        "MRE11A"        "MTX2"          "HADHA"         "NRP2"         
## [185] "LYPLA2"        "SHKBP1"        "GLOD4"         "ATP6V1H"       "EIF4E2"        "NME1"          "HDAC2"         "JUP"          
## [193] "PHLDB1"        "NUTF2"         "CCT7"          "STRAP"         "TES"           "BUB3"          "DLX5"          "SYF2"         
## [201] "GPBP1L1"       "DDB1"          "EZH2"          "SNRPA1"        "PDCD6"         "EMC2"          "RP4-702J19.1"  "CAV1"         
## [209] "CDK4"          "CCT5"          "RAB1A"         "CDKN1A"        "ELAVL1"        "MBD4"          "PLEKHJ1"       "CCT3"         
## [217] "ATPIF1"        "CDCA7"         "C22orf29"      "SEC23B"        "SPIDR"         "HSPB8"         "UBE2I"         "CPNE3"        
## [225] "NEFH"          "LMAN1"         "RAB18"         "ILK"           "EIF2A"         "SLBP"          "ELOVL1"        "MYB"          
## [233] "RP4-753D5.4"   "KRT18"         "CRTAP"         "MUM1"          "TPI1"          "KIAA1715"      "RP11-259K5.3"  "DTYMK"        
## [241] "MAEA"          "PNRC2"         "GLT25D1"       "MCM6"          "RBM17"         "TTC39A"        "FEN1"          "SHARPIN"      
## [249] "IFITM1"        "PSMD4"         "SNORD99"       "BZW2"          "TIMMDC1"       "STAG3"         "TPD52"         "SNRPG"        
## [257] "TIMM13"        "JAG1"          "BTG2"          "LAMTOR1"       "AC019188.1"    "VAMP3"         "CAD"           "ATAD2"

myb_targets <- read_excel(input$myb_targets,col_names = T) %>%  pull(Gene)

tp53_data = read_tsv(file = input$tp53_data,col_names = T)

## Rows: 530 Columns: 8
## [36mℹ[39m Targets are out of date. Updating...

## [36mℹ[39m Targets are out of date. Updating...── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## [36mℹ[39m Targets are out of date. Updating...Delimiter: "\t"
## chr (7): Study ID, Sample ID, Patient ID, TP53, TP53: MUT, TP53: AMP, TP53: HOMDEL
## dbl (1): Altered
## 
## [36mℹ[39m Targets are out of date. Updating...
## 
## [36mℹ[39m Targets are out of date. Updating...ℹ Use `spec()` to retrieve the full column specification for this data.
## [36mℹ[39m Targets are out of date. Updating...ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## [36mℹ[39m Targets are out of date. Updating...

# load TPM data
HNSC_tpm = readRDS(input$HNSC_tpm)

if (params$log_tpm) {
  HNSC_tpm = log2(HNSC_tpm+1)
}

# load clinical data
project_name = "TCGA-HNSC"

clinical_data <- GDCquery_clinic(project_name, "clinical")
OPSCC_tissues = c("Base of tongue, NOS", #filter only OPSCC
"Oropharynx, NOS",
"Posterior wall of oropharynx",
"Tonsil, NOS")
clinical_data = clinical_data[clinical_data$tissue_or_organ_of_origin %in% OPSCC_tissues,]

#add follow_up data
clinical_follow_up = read_tsv(input$clinical_follow_up) %>% as.data.frame()

## New names:
## Rows: 1949 Columns: 139 [36mℹ[39m Targets are out of date. Updating...
## [36mℹ[39m Targets are out of date. Updating...── Column specification
## ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── [36mℹ[39m Targets are out of
## date. Updating...Delimiter: "\t" chr (139): case_id, case_submitter_id, project_id, adverse_event, adverse_event_grade, aids_risk_factors,
## barretts_esophagus_goblet_cells_prese...
## 
## [36mℹ[39m Targets are out of date. Updating...
## [36mℹ[39m Targets are out of date. Updating...ℹ Use `spec()` to retrieve the full column specification for this data. [36mℹ[39m Targets are out of date.
## Updating...ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message. [36mℹ[39m Targets are out of date. Updating...
## • `timepoint_category` -> `timepoint_category...79`
## • `timepoint_category` -> `timepoint_category...135`

#get only "Last Contact"
clinical_follow_up <- clinical_follow_up[clinical_follow_up$timepoint_category...135 == "Last Contact", ]

#add to clinical and rename days_to_last_follow_up> days_to_follow_up
clinical_data = inner_join(clinical_data,clinical_follow_up[,c("case_submitter_id","days_to_follow_up"),drop=F],by = c("submitter_id"="case_submitter_id")) %>% select(-days_to_last_follow_up) %>% dplyr::rename(days_to_last_follow_up = days_to_follow_up)

# add HPV statusfrom cbioportal to clinical_data
cbp_data = read_tsv(file = input$cbp_data)

## Rows: 523 Columns: 62
## [36mℹ[39m Targets are out of date. Updating...
## [36mℹ[39m Targets are out of date. Updating...── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## [36mℹ[39m Targets are out of date. Updating...Delimiter: "\t"
## chr (43): Study ID, Patient ID, Sample ID, Neoplasm Disease Stage American Joint Committee on Cancer Code, American Joint Committee on Cancer P...
## dbl (18): Diagnosis Age, Aneuploidy Score, Buffa Hypoxia Score, Last Communication Contact from Initial Pathologic Diagnosis Date, Birth from I...
## lgl  (1): Patient Weight
## 
## [36mℹ[39m Targets are out of date. Updating...
## 
## [36mℹ[39m Targets are out of date. Updating...ℹ Use `spec()` to retrieve the full column specification for this data.
## [36mℹ[39m Targets are out of date. Updating...ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## [36mℹ[39m Targets are out of date. Updating...

cbp_data = cbp_data[,c("Patient ID", "Subtype")] %>% dplyr::rename(hpv_status = Subtype) %>% 
  mutate(hpv_status = gsub(,x = hpv_status,pattern = "HNSC_HPV",replacement = "HPV"))
clinical_data = inner_join(x = clinical_data,y = cbp_data,by = c("submitter_id"="Patient ID"))
na_patients = clinical_data[is.na(clinical_data$hpv_status), "submitter_id"]

# add TP53 mutation status from cbioportal to clinical_data

tp53_data = tp53_data[,c("Patient ID", "Altered")] %>% dplyr::rename(tp53_status = Altered) %>% 
  mutate(tp53_status = 
           case_when(tp53_status == 1 ~ "tp53_mutated",
                     tp53_status == 0 ~ "tp53_unmutated"))
clinical_data = inner_join(x = clinical_data,y = tp53_data,by = c("submitter_id"="Patient ID"))
na_patients = clinical_data[is.na(clinical_data$tp53_data), "submitter_id"]

# #set clinical data to max 5 years
# clinical_data = clinical_data %>% mutate(days_to_last_follow_up = as.numeric(days_to_last_follow_up)) %>% 
#   mutate(vital_status = if_else(condition = (days_to_death) > 1825 | is.na(days_to_death) , false = vital_status, true = "Alive")) %>% 
#   mutate(days_to_last_follow_up = if_else(condition = days_to_last_follow_up > 1825, false = days_to_last_follow_up,true = 1825)) %>%
#   mutate(days_to_death = if_else(condition = days_to_death > 1825, false = days_to_death, true = NA))

2 MYB

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
paged_table(gene_status)

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))


p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p[[1]]/p[[2]]

pdf(paste0(params$data_out_dir,"KM_MYB_OPSCC.pdf"),width = 8,height = 6)
p[[1]]/p[[2]]
dev.off()

## png 
##   2

3 HPV

clinical_data_for_hpv = clinical_data  %>% filter(!is.na(hpv_status))
p = TCGAanalyze_survival(
    data = clinical_data_for_hpv,
    clusterCol = "hpv_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p[[1]]/p[[2]]

4 MYB - HPV PPM correlation

# load data 
library("readxl")
library(stringr)
library(ggpubr)
TCGA_HPV_data <- read_excel("./Input_data/TCGA/PMC3806554_ncomms3513-s5.xlsx",skip = 2)
TCGA_HPV_data$submitter_id = TCGA_HPV_data$`Sample Barcode` %>% str_sub( start = 1, end = 12) 
rownames(TCGA_HPV_data) = TCGA_HPV_data$submitter_id

## Warning: Setting row names on a tibble is deprecated.

TCGA_HPV_data = TCGA_HPV_data %>% filter(Cancer == "HNSC")

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status = gene_status[TCGA_HPV_data$submitter_id,] %>%  cbind(HPV_PPM = TCGA_HPV_data$ppm) %>% na.omit()

paged_table(gene_status)

sp <- ggscatter(gene_status, x = "MYB", y = "HPV_PPM",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   )
# Add correlation coefficient
sp + stat_cor(method = "pearson")

## `geom_smooth()` using formula = 'y ~ x'

5 HPV signature - HPV correlation

# with HPV signature
gene_status = set_clinical_data(clin_data = clinical_data,genes = opscc_signature,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HPV+ signature")
gene_status = gene_status[TCGA_HPV_data$submitter_id,] %>% cbind(HPV_PPM = TCGA_HPV_data$ppm) %>% na.omit()
paged_table(gene_status)

sp <- ggscatter(gene_status, x = "HPV+ signature", y = "HPV_PPM",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   )
# Add correlation coefficient
sp + stat_cor(method = "pearson")

## `geom_smooth()` using formula = 'y ~ x'

6 HPV signature - HPV T test

# create data

gene_status = set_clinical_data(clin_data = clinical_data,genes = opscc_signature,tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)
gene_status = inner_join(x = gene_status,y = clinical_data[,c("submitter_id","hpv_status")],by = c("patient"="submitter_id"))

paged_table(gene_status)

# plot
library(ggpubr)
library(rstatix)
stat.test <- gene_status %>%
  t_test(score ~ hpv_status) %>%
  add_significance()
stat.test

bxp <- ggboxplot(gene_status, x = "hpv_status", y = "score", fill = "#00AFBB")+
  geom_jitter()
stat.test <- stat.test %>% add_xy_position(x = "hpv_status")
bxp + 
  stat_pvalue_manual(stat.test,  label = "T-test, p = {p}")

7 MYB- HPV T test

# create data

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)
gene_status = inner_join(x = gene_status,y = clinical_data[,c("submitter_id","hpv_status")],by = c("patient"="submitter_id"))

paged_table(gene_status)

# plot
library(ggpubr)
library(rstatix)
gene_status$hpv_status = factor(gene_status$hpv_status,levels = c("HPV-","HPV+"))
stat.test <- gene_status %>%
  t_test(MYB ~ hpv_status) %>%
  add_significance()
stat.test

stat.test <- stat.test %>% add_xy_position(x = "hpv_status")


plt = ggplot(gene_status, aes(x = hpv_status, y = MYB)) +
  geom_violin(scale = "width",aes(fill = hpv_status)) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_pvalue_manual(stat.test, label = "{p.adj}")  + #add p value
  geom_jitter()+
  labs(title = "MYB expression by HPV status", y = "MYB expression: log2(TPM)", fill = "HPV Status")+
  theme(plot.title = element_text(hjust = 0.5),axis.title.x = element_blank())+scale_x_discrete(limits = c("HPV-","HPV+"))+theme_minimal()
plt

## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_ydensity()`).

## Warning: Removed 5 rows containing missing values or values outside the scale range (`stat_boxplot()`).

## Warning: Removed 5 rows containing missing values or values outside the scale range (`geom_point()`).

pdf(paste0(params$data_out_dir,"violin_MYB_HPV_OPSCC_TCGA.pdf"),width = 7,height = 5)
plt

## Warning: Removed 5 rows containing non-finite outside the scale range (`stat_ydensity()`).

## Warning: Removed 5 rows containing missing values or values outside the scale range (`stat_boxplot()`).

## Warning: Removed 5 rows containing missing values or values outside the scale range (`geom_point()`).

dev.off()

## png 
##   2

8 MYB and HPV survival plot

8.1 HPV + tumors

clinical_data_with_scores = clinical_data  %>% filter(hpv_status == "HPV+")

gene_status = set_clinical_data(clin_data = clinical_data_with_scores,genes = "MYB",tpm_data_frame = HNSC_tpm,stratify = "T")
paged_table(gene_status)

gene_status$patient = rownames(gene_status)
clinical_data_with_scores = inner_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)
p[[1]]/p[[2]]

8.2 HPV- tumors

clinical_data_with_scores = clinical_data  %>% filter(hpv_status == "HPV-")

gene_status = set_clinical_data(clin_data = clinical_data_with_scores,genes = "MYB",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
paged_table(gene_status)

gene_status$patient = rownames(gene_status)
clinical_data_with_scores = inner_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)
p[[1]]/p[[2]]

8.3 combine MYB and HPV survival plot

gene_status = set_clinical_data(clin_data = clinical_data,genes = "MYB",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)

clinical_data_with_scores = inner_join(x = clinical_data,y = gene_status, by = c("submitter_id"="patient"))

clinical_data_with_scores$gene_status = paste(clinical_data_with_scores$gene_status,clinical_data_with_scores$hpv_status,sep = ", ")
paged_table(clinical_data_with_scores)

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)
p[[1]]/p[[2]]

9 Myb Targets

# create data

gene_status = set_clinical_data(clin_data = clinical_data,genes = "NOTCH1",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)
gene_status = inner_join(x = gene_status,y = clinical_data[,c("submitter_id","hpv_status")],by = c("patient"="submitter_id"))

paged_table(gene_status)

# plot
library(ggpubr)
library(rstatix)
stat.test <- gene_status %>%
  t_test(NOTCH1 ~ hpv_status) %>%
  add_significance()
stat.test

bxp <- ggboxplot(gene_status, x = "hpv_status", y = "NOTCH1", fill = "#00AFBB")+
  geom_jitter()
stat.test <- stat.test %>% add_xy_position(x = "hpv_status")
bxp + 
  stat_pvalue_manual(stat.test,  label = "T-test, p = {p}") +ylim(0,7)

10 HMSC signature

# all samples:
gene_status = set_clinical_data(clin_data = clinical_data,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HPV+ signature")

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p1 = (p[[1]]+ggtitle("OPSCC samples"))/p[[2]]
print_tab(p1,title = "all samples")

all samples

pdf(paste0(params$data_out_dir,"all_samples_hmsc_signature.pdf"),width = 8,height = 6)
p1
dev.off()

png 2

# HPV+ samples:
clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV+")
gene_status = set_clinical_data(clin_data = clinical_data_hpv_pos,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature")
paged_table(gene_status)

clinical_data_with_scores = clinical_data_hpv_pos[clinical_data_hpv_pos$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p1 = (p[[1]]+ggtitle("OPSCC HPV+ samples"))/p[[2]]
print_tab(p1,title = "HPV+")

HPV+

pdf(paste0(params$data_out_dir,"HPVpos_samples_hmsc_signature.pdf"),width = 8,height = 6)
p1
dev.off()

png 2

# HPV neg samples:
clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV-")
gene_status = set_clinical_data(clin_data = clinical_data_hpv_pos,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature")
clinical_data_with_scores = clinical_data_hpv_pos[clinical_data_hpv_pos$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p1 = (p[[1]]+ggtitle("OPSCC HPV- samples"))/p[[2]]
print_tab(p1,title = "HPV-")

HPV-

pdf(paste0(params$data_out_dir,"HPVneg_samples_hmsc_signature.pdf"),width = 8,height = 6)
p1
dev.off()

png 2

10.1 No cell cycle genes

# remove genes:
genesets_h  = geneIds(getGmt(input$genesets_h))
message("cc genes in signature:")

## cc genes in signature:

print(hmsc_signatue[hmsc_signatue %in% genesets_h$GO_MITOTIC_CELL_CYCLE])

[1] “ANLN” “TUBB6” “PBRM1” “MCM5” “TUBA1C” “MSH6” “KNTC1” “SMC5” “DCTN2” “IK” “TAOK3” “PSME2” “SEH1L” “MCM3” “PSMB1” [16] “MCM4” “PCNT” “KIF22” “PSMB4” “VPS4B” “CNOT7” “SMC3” “TYMS” “PSMA5” “MYC” “BUB3” “SYF2” “DDB1” “EZH2” “CDK4”
[31] “CDKN1A” “MYB” “MAEA” “MCM6” “PSMD4” “BTG2”

hmsc_signatue_no_cc = hmsc_signatue[!hmsc_signatue %in% genesets_h$GO_MITOTIC_CELL_CYCLE]

#all samples:

gene_status = set_clinical_data(clin_data = clinical_data,genes = hmsc_signatue_no_cc,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature no CC")

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p1 = (p[[1]]+ggtitle("OPSCC samples"))/p[[2]]
print_tab(p1,title = "OPSCC samples",subtitle_num = 3)

OPSCC samples

pdf(paste0(params$data_out_dir,"all_samples_hmsc_signature_noCC_genes.pdf"),width = 8,height = 6)
p1
dev.off()

png 2

#hpv +:
clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV+")

gene_status = set_clinical_data(clin_data = clinical_data_hpv_pos,genes = hmsc_signatue_no_cc,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature no CC")

clinical_data_with_scores = clinical_data_hpv_pos[clinical_data_hpv_pos$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p1 = (p[[1]]+ggtitle("OPSCC HPV+ samples"))/p[[2]]
print_tab(p1,title = "HPV+",subtitle_num = 3)

HPV+

pdf(paste0(params$data_out_dir,"HPVpos_samples_hmsc_signature_noCC_genes.pdf"),width = 8,height = 6)
p1
dev.off()

png 2

#hpv- :
clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV-")

gene_status = set_clinical_data(clin_data = clinical_data_hpv_pos,genes = hmsc_signatue_no_cc,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature no CC")

clinical_data_with_scores = clinical_data_hpv_pos[clinical_data_hpv_pos$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p1 = (p[[1]]+ggtitle("OPSCC HPV- samples"))/p[[2]]
print_tab(p1,title = "HPV-",subtitle_num = 3)

HPV-

pdf(paste0(params$data_out_dir,"HPVneg_samples_hmsc_signature_noCC_genes.pdf"),width = 8,height = 6)
p1
dev.off()

png 2

10.2 T test

# create data

gene_status = set_clinical_data(clin_data = clinical_data,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)
gene_status = inner_join(x = gene_status,y = clinical_data[,c("submitter_id","hpv_status")],by = c("patient"="submitter_id"))
gene_status = gene_status %>% filter(!is.na(hpv_status))
paged_table(gene_status)

# plot
library(ggpubr)
library(rstatix)
stat.test <- gene_status %>%
  t_test(score ~ hpv_status) %>%
  add_significance()
stat.test

bxp <- ggboxplot(gene_status, x = "hpv_status", y = "score", fill = "#00AFBB")+
  geom_jitter()
stat.test <- stat.test %>% add_xy_position(x = "hpv_status")
bxp + 
  stat_pvalue_manual(stat.test,  label = "T-test, p = {p}")

# HPV signature ~ stage

gene_status = set_clinical_data(clin_data = clinical_data,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify)
gene_status$patient = rownames(gene_status)
gene_status = inner_join(x = gene_status,y = clinical_data[,c("submitter_id","ajcc_clinical_stage","hpv_status")],by = c("patient"="submitter_id"))
gene_status = gene_status %>% filter(!is.na(ajcc_clinical_stage))
paged_table(gene_status)

gene_status$ajcc_clinical_stage = factor(gene_status$ajcc_clinical_stage,levels = c("Stage I", "Stage II", "Stage III","Stage IVA", "Stage IVB", 
"Stage IVC"))

my_comparisons <- list( c("Stage I", "Stage II"),c("Stage II", "Stage III"),c("Stage III", "Stage IVA"),c("Stage IVA", "Stage IVB"),c("Stage IVB", "Stage IVC") )

p_list = list()
for (subset_name in list("HPV+","HPV-",c("HPV+","HPV-"))) {
  gene_status_subset = gene_status %>% filter(hpv_status %in% subset_name)
  
  stat.test <-  
  compare_means(score ~ ajcc_clinical_stage,data = gene_status_subset,
                method = "wilcox.test") %>%
  filter(paste(group1, group2) %in% unlist(lapply(my_comparisons, paste, collapse = " "))) %>% 
    mutate(p = round(p, 3))

  stat.test = stat.test %>% add_y_position(formula = score ~ ajcc_clinical_stage,data = gene_status_subset,comparisons = my_comparisons,step.increase = 0)
  
  p <- ggboxplot(gene_status_subset, x = "ajcc_clinical_stage", y = "score",
                 color = "ajcc_clinical_stage", palette = "npg",
                add = "jitter",)+
      stat_pvalue_manual(stat.test, label = "p = {p}",bracket.shorten = 0.1)+
      ylab("HMSC_HPV_score")+
    ggtitle("OPSCC TCGA " %s+% paste0(subset_name,collapse = " and ") %s+% " samples")+
    theme(legend.position = "none")
  p_list[[paste0(subset_name,collapse = " and ")]] = p
  print(p)
}

pdf(file = params$data_out_dir %s+% "HPV_signature_score-stage_OPSCC_TCGA.pdf")
lapply(p_list,function(x) {x})

## $`HPV+`

## 
## $`HPV-`

## 
## $`HPV+ and HPV-`

dev.off()

## png 
##   2

11 E2F ~ MYB

clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV+")
clinical_data_hpv_neg = clinical_data  %>% filter(hpv_status == "HPV-")

all_figs = list()
for (gene in c("E2F1","E2F2","E2F3")) {
  opscc_tpm = HNSC_tpm[,colnames(HNSC_tpm) %in% clinical_data_hpv_pos$submitter_id] #filter for OPSCC tissues as in clinical data
  opscc_tpm = opscc_tpm[c("MYB",gene),] %>% t()
  
  sp1 <- ggscatter(opscc_tpm, x = "MYB", y = gene,
     add = "reg.line",  # Add regressin line
     add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
     conf.int = TRUE # Add confidence interval
     ) +
    stat_cor(method = "pearson")+ # Add correlation coefficient
    ggtitle("OPSCC TCGA HPV+ samples")+
    ylab(gene %s+% " log(TPM)")+
    xlab("MYB log(TPM)")
  sp1
  
  
  # with seperate HPV- cell analysis
  # opscc_tpm = HNSC_tpm[,colnames(HNSC_tpm) %in% clinical_data_hpv_neg$submitter_id] #filter for OPSCC tissues as in clinical data
  # opscc_tpm = opscc_tpm[c("MYB",gene),] %>% t()
  # 
  # sp2 <- ggscatter(opscc_tpm, x = "MYB", y = gene,
  #    add = "reg.line",  # Add regressin line
  #    add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
  #    conf.int = TRUE # Add confidence interval
  #    ) +
  #   stat_cor(method = "pearson")+ # Add correlation coefficient
  #   ggtitle("OPSCC TCGA HPV- samples")+
  #   ylab(gene %s+% " log(TPM)")+
  #   xlab("MYB log(TPM)")
  # sp2
  # all_figs[[gene]] = sp1/sp2
  all_figs[[gene]] = sp1

}

pdf(file = params$data_out_dir %s+% "TCGA MYB-ETF correlation_HPV_pos.pdf",height = 5,width = 6)
lapply(all_figs,function(x) {x})

## $E2F1

## `geom_smooth()` using formula = 'y ~ x'

## 
## $E2F2

## `geom_smooth()` using formula = 'y ~ x'

## 
## $E2F3

## `geom_smooth()` using formula = 'y ~ x'

dev.off()

## png 
##   2

12 MYB ~ RBL1

clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV+")

opscc_tpm = HNSC_tpm[,colnames(HNSC_tpm) %in% clinical_data_hpv_pos$submitter_id] #filter for OPSCC tissues as in clinical data
opscc_tpm = opscc_tpm[c("MYB","RBL1"),] %>% t()

sp <- ggscatter(opscc_tpm, x = "MYB", y = "RBL1",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   ) +
  stat_cor(method = "pearson")+ # Add correlation coefficient
  ggtitle("OPSCC TCGA samples")+
  ylab("E2F1 log(TPM)")+
  xlab("MYB log(TPM)")
sp

## `geom_smooth()` using formula = 'y ~ x'

clinical_data_hpv_neg = clinical_data  %>% filter(hpv_status == "HPV-")

opscc_tpm = HNSC_tpm[,colnames(HNSC_tpm) %in% clinical_data_hpv_neg$submitter_id] #filter for OPSCC tissues as in clinical data
opscc_tpm = opscc_tpm[c("MYB","RBL1"),] %>% t()

sp <- ggscatter(opscc_tpm, x = "MYB", y = "RBL1",
   add = "reg.line",  # Add regressin line
   add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
   conf.int = TRUE # Add confidence interval
   ) +
  stat_cor(method = "pearson")+ # Add correlation coefficient
  ggtitle("OPSCC TCGA samples")+
  ylab("E2F1 log(TPM)")+
  xlab("MYB log(TPM)")
sp

## `geom_smooth()` using formula = 'y ~ x'

13 OPSCC sc signature

13.1 all samples

gene_status = set_clinical_data(clin_data = clinical_data,genes = opscc_signature,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "OPSCC HPV \nsignature")
paged_table(gene_status)

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

(p[[1]]+ggtitle("OPSCC samples"))/p[[2]]

13.2 HPV +

clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV+")

gene_status = set_clinical_data(clin_data = clinical_data_hpv_pos,genes = opscc_signature,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "OPSCC HPV \nsignature")
paged_table(gene_status)

clinical_data_with_scores = clinical_data_hpv_pos[clinical_data_hpv_pos$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

(p[[1]]+ggtitle("OPSCC HPV+ samples"))/p[[2]]

13.3 HPV -

clinical_data_hpv_pos = clinical_data  %>% filter(hpv_status == "HPV-")

gene_status = set_clinical_data(clin_data = clinical_data_hpv_pos,genes = opscc_signature,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "OPSCC HPV \nsignature")
paged_table(gene_status)

clinical_data_with_scores = clinical_data_hpv_pos[clinical_data_hpv_pos$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

(p[[1]]+ggtitle("OPSCC HPV- samples"))/p[[2]]

14 Additional mutations

gene_status = set_clinical_data(clin_data = clinical_data,genes = "TP53",tpm_data_frame = HNSC_tpm,stratify = params$stratify)
paged_table(gene_status)

clinical_data_with_scores = clinical_data[clinical_data$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    main = "TCGA Set\n GBM",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p[[1]]/p[[2]]

14.1 TP53

p = TCGAanalyze_survival(
    data = clinical_data,
    clusterCol = "tp53_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)

p[[1]]/p[[2]]

14.2 Fisher

# plot
library(rstatix)

df = clinical_data[,c("tp53_status","hpv_status")]
test = fisher_test(table(df))
p = ggbarstats(
  df,
  tp53_status,
  hpv_status,
  results.subtitle = FALSE,
  subtitle = paste0("Fisher's exact test", ", p-value = ",
                    test$p)
)

p

14.3 TP53 unmutated

clinical_data_for_test = clinical_data  %>% filter(tp53_status == "tp53_unmutated")

gene_status = set_clinical_data(clin_data = clinical_data_for_test,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature")
paged_table(gene_status)

clinical_data_with_scores = clinical_data_for_test[clinical_data_for_test$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)
(p[[1]]+ggtitle("OPSCC unmutated tp53 samples"))/p[[2]]

14.4 TP53 mutated

clinical_data_for_test = clinical_data  %>% filter(tp53_status == "tp53_mutated")

gene_status = set_clinical_data(clin_data = clinical_data_for_test,genes = hmsc_signatue,tpm_data_frame = HNSC_tpm,stratify = params$stratify,signature_name = "HMSC HPV \nsignature") 
paged_table(gene_status)

clinical_data_with_scores = clinical_data_for_test[clinical_data_for_test$submitter_id %in% rownames(gene_status),]
gene_status$patient = rownames(gene_status)
clinical_data_with_scores = left_join(x = clinical_data_with_scores,y = gene_status, by = c("submitter_id"="patient"))

p = TCGAanalyze_survival(
    data = clinical_data_with_scores,
    clusterCol = "gene_status",
    height = 10, width=10,filename = NULL,
    pval.method =T ,
    xscale = "d_m",break.time.by = 365.2422,
    xlab = "Time since diagnosis (months)"
)
(p[[1]]+ggtitle("OPSCC mutated tp53 samples"))/p[[2]]

15 session info

library(devtools)
session_info()

## ─ Session info ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.1 (2024-06-14)
##  os       Ubuntu 22.04.5 LTS
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Etc/UTC
##  date     2025-09-02
##  pandoc   3.4 @ /usr/bin/ (via rmarkdown)
## 
## ─ Packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  !  package              * version    date (UTC) lib source
##     abind                  1.4-5      2016-07-21 [1] RSPM (R 4.1.0)
##     annotate             * 1.84.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     AnnotationDbi        * 1.56.2     2021-11-09 [1] Bioconductor
##     ape                    5.8        2024-04-11 [1] RSPM (R 4.4.0)
##     aplot                  0.1.3      2022-04-01 [1] RSPM (R 4.1.0)
##     babelgene              22.3       2022-03-30 [1] RSPM (R 4.1.0)
##     backports              1.5.0      2024-05-23 [1] RSPM (R 4.4.0)
##     bayestestR             0.15.0     2024-10-17 [1] RSPM (R 4.4.0)
##     Biobase              * 2.54.0     2021-10-26 [1] Bioconductor
##     BiocFileCache          2.2.1      2022-01-23 [1] Bioconductor
##     BiocGenerics         * 0.52.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     BiocParallel           1.40.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     biomaRt              * 2.50.3     2022-02-03 [1] Bioconductor
##     Biostrings             2.62.0     2021-10-26 [1] Bioconductor
##     bit                    4.0.4      2020-08-04 [1] RSPM
##     bit64                  4.0.5      2020-08-30 [1] RSPM
##     bitops                 1.0-7      2021-04-24 [1] RSPM
##     blob                   1.2.3      2022-04-10 [1] RSPM
##     broom                  1.0.7      2024-09-26 [1] RSPM (R 4.4.0)
##     bslib                  0.9.0      2025-01-30 [1] RSPM (R 4.4.0)
##     cachem                 1.0.6      2021-08-19 [1] RSPM
##     car                    3.0-12     2021-11-06 [1] RSPM (R 4.1.0)
##     carData                3.0-5      2022-01-06 [1] RSPM (R 4.4.0)
##     cellranger             1.1.0      2016-07-27 [1] RSPM (R 4.4.0)
##     chromote               0.3.1      2024-08-30 [1] RSPM (R 4.4.0)
##     circlize               0.4.16     2024-02-20 [1] RSPM (R 4.4.0)
##     cli                    3.6.3      2024-06-21 [1] RSPM (R 4.4.0)
##     clue                   0.3-66     2024-11-13 [1] RSPM (R 4.4.0)
##     cluster                2.1.6      2023-12-01 [3] CRAN (R 4.4.1)
##     clusterProfiler      * 4.2.2      2022-01-13 [1] Bioconductor
##     codetools              0.2-20     2024-03-31 [3] CRAN (R 4.4.1)
##     colorspace             2.0-3      2022-02-21 [1] RSPM
##     ComplexHeatmap       * 2.10.0     2021-10-26 [1] Bioconductor
##     conflicted           * 1.2.0      2023-02-01 [1] RSPM (R 4.4.0)
##     correlation            0.8.6      2024-10-26 [1] RSPM (R 4.4.0)
##     cowplot                1.1.1      2020-12-30 [1] RSPM (R 4.1.0)
##     crayon                 1.5.1      2022-03-26 [1] RSPM
##     curl                   6.0.0      2024-11-05 [1] RSPM (R 4.4.0)
##     data.table           * 1.14.2     2021-09-27 [1] RSPM (R 4.1.0)
##     datawizard             0.13.0     2024-10-05 [1] RSPM (R 4.4.0)
##     DBI                    1.2.3      2024-06-02 [1] RSPM (R 4.4.0)
##     dbplyr                 2.3.4      2023-09-26 [1] RSPM (R 4.4.0)
##     DelayedArray           0.24.0     2022-11-01 [1] Bioconductor
##     deldir                 1.0-6      2021-10-23 [1] RSPM (R 4.1.0)
##     devtools             * 2.4.3      2021-11-30 [1] RSPM (R 4.1.0)
##     digest                 0.6.29     2021-12-01 [1] RSPM
##     DO.db                  2.9        2022-11-08 [1] Bioconductor
##     doParallel             1.0.17     2022-02-07 [1] RSPM
##     DOSE                   3.20.1     2021-11-18 [1] Bioconductor
##     downloader             0.4        2015-07-09 [1] RSPM (R 4.1.0)
##     dplyr                * 1.1.4      2023-11-17 [1] RSPM (R 4.4.0)
##     effectsize             0.8.9      2024-07-03 [1] RSPM (R 4.4.0)
##     ellipsis               0.3.2      2021-04-29 [1] RSPM
##     enrichplot             1.14.2     2022-02-24 [1] Bioconductor
##     evaluate               0.15       2022-02-18 [1] RSPM
##     fansi                  1.0.3      2022-03-24 [1] RSPM
##     farver                 2.1.0      2021-02-28 [1] RSPM
##     fastmap                1.1.1      2023-02-24 [1] RSPM (R 4.4.0)
##     fastmatch              1.1-4      2023-08-18 [1] RSPM (R 4.4.0)
##     fdrtool              * 1.2.17     2021-11-13 [1] RSPM (R 4.1.0)
##     fgsea                  1.20.0     2021-10-26 [1] Bioconductor
##     filelock               1.0.3      2023-12-11 [1] CRAN (R 4.2.2)
##     fitdistrplus           1.1-8      2022-03-10 [1] RSPM (R 4.1.0)
##     forcats              * 1.0.0      2023-01-29 [1] RSPM (R 4.4.0)
##     foreach                1.5.2      2022-02-02 [1] RSPM
##     fs                     1.5.2      2021-12-08 [1] RSPM
##     future                 1.24.0     2022-02-19 [1] RSPM (R 4.1.0)
##     future.apply           1.8.1      2021-08-10 [1] RSPM (R 4.1.0)
##     gargle                 1.5.2      2023-07-20 [1] RSPM (R 4.4.0)
##     generics               0.1.2      2022-01-31 [1] RSPM
##     GenomeInfoDb           1.30.1     2022-01-30 [1] Bioconductor
##     GenomeInfoDbData       1.2.7      2024-01-17 [1] Bioconductor
##     GenomicRanges          1.58.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     GetoptLong             1.0.5      2020-12-15 [1] RSPM (R 4.4.0)
##     ggforce                0.3.3      2021-03-05 [1] RSPM
##     ggfun                  0.0.6      2022-04-01 [1] RSPM (R 4.1.0)
##     ggplot2              * 3.5.1      2024-04-23 [1] RSPM (R 4.4.0)
##     ggplotify              0.1.0      2021-09-02 [1] RSPM (R 4.1.0)
##     ggpubr               * 0.4.0      2020-06-27 [1] RSPM (R 4.1.0)
##     ggraph                 2.0.5      2021-02-23 [1] RSPM (R 4.1.0)
##     ggrepel              * 0.9.6      2024-09-07 [1] RSPM (R 4.4.0)
##     ggridges               0.5.3      2021-01-08 [1] RSPM (R 4.1.0)
##     ggsci                  2.9        2018-05-14 [1] RSPM (R 4.1.0)
##     ggsignif               0.6.4      2022-10-13 [1] RSPM (R 4.4.0)
##     ggstatsplot          * 0.12.5     2024-11-01 [1] RSPM (R 4.4.0)
##     ggtext                 0.1.1      2020-12-17 [1] RSPM (R 4.1.0)
##     ggtree                 3.2.1      2021-11-16 [1] Bioconductor
##     GlobalOptions          0.1.2      2020-06-10 [1] RSPM (R 4.4.0)
##     globals                0.14.0     2020-11-22 [1] RSPM (R 4.1.0)
##     glue                   1.8.0      2024-09-30 [1] RSPM (R 4.4.0)
##     GO.db                  3.14.0     2022-11-08 [1] Bioconductor
##     goftest                1.2-3      2021-10-07 [1] RSPM (R 4.1.0)
##     googledrive            2.1.1      2023-06-11 [1] RSPM (R 4.4.0)
##     googlesheets4          1.1.1      2023-06-11 [1] RSPM (R 4.4.0)
##     GOSemSim               2.20.0     2021-10-26 [1] Bioconductor
##     graph                * 1.72.0     2021-10-26 [1] Bioconductor
##     graphlayouts           0.8.0      2022-01-03 [1] RSPM (R 4.1.0)
##     gridExtra              2.3        2017-09-09 [1] RSPM
##     gridGraphics           0.5-1      2020-12-13 [1] RSPM (R 4.1.0)
##     gridtext               0.1.4      2020-12-10 [1] RSPM (R 4.1.0)
##     GSEABase             * 1.56.0     2021-10-26 [1] Bioconductor
##     gtable                 0.3.0      2019-03-25 [1] RSPM
##     haven                  2.5.4      2023-11-30 [1] RSPM (R 4.4.0)
##     hms                    1.1.3      2023-03-21 [1] RSPM (R 4.4.0)
##     htmltools              0.5.8.1    2024-04-04 [1] RSPM (R 4.4.0)
##     htmlwidgets            1.5.4      2021-09-08 [1] RSPM (R 4.1.0)
##     httpuv                 1.6.5      2022-01-05 [1] RSPM (R 4.1.0)
##     httr                   1.4.7      2023-08-15 [1] RSPM (R 4.4.0)
##     ica                    1.0-2      2018-05-24 [1] RSPM (R 4.1.0)
##  VP igraph               * 1.4.2      2024-10-19 [?] RSPM (R 4.4.0) (on disk 2.1.1)
##     insight                0.20.5     2024-10-02 [1] RSPM (R 4.4.0)
##     IRanges              * 2.40.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     irlba                  2.3.5      2021-12-06 [1] RSPM (R 4.1.0)
##     iterators              1.0.14     2022-02-05 [1] RSPM
##     job                    0.3.1      2024-05-05 [1] RSPM (R 4.4.0)
##     jquerylib              0.1.4      2021-04-26 [1] RSPM
##     jsonlite               1.8.9      2024-09-20 [1] RSPM (R 4.4.0)
##     kableExtra           * 1.3.4      2021-02-20 [1] RSPM (R 4.1.0)
##     KEGGREST               1.34.0     2021-10-26 [1] Bioconductor
##     KernSmooth             2.23-24    2024-05-17 [3] CRAN (R 4.4.1)
##     km.ci                  0.5-6      2022-04-06 [1] RSPM (R 4.1.0)
##     KMsurv                 0.1-5      2012-12-03 [1] RSPM (R 4.1.0)
##     knitr                  1.49       2024-11-08 [1] RSPM (R 4.4.0)
##     labeling               0.4.2      2020-10-20 [1] RSPM
##     later                  1.3.0      2021-08-18 [1] RSPM (R 4.1.0)
##     lattice                0.22-6     2024-03-20 [3] CRAN (R 4.4.1)
##     lazyeval               0.2.2      2019-03-15 [1] RSPM (R 4.1.0)
##     leiden                 0.3.9      2021-07-27 [1] RSPM (R 4.1.0)
##     lifecycle              1.0.4      2023-11-07 [1] RSPM (R 4.4.0)
##     listenv                0.8.0      2019-12-05 [1] RSPM (R 4.1.0)
##     lmtest                 0.9-40     2022-03-21 [1] RSPM (R 4.1.0)
##     lubridate              1.9.3      2023-09-27 [1] RSPM (R 4.4.0)
##     magrittr             * 2.0.3      2022-03-30 [1] RSPM (R 4.4.0)
##     makepipe             * 0.2.2      2025-01-07 [1] RSPM (R 4.4.0)
##     markdown               1.1        2019-08-07 [1] RSPM (R 4.1.0)
##     MASS                   7.3-60.2   2024-04-26 [3] CRAN (R 4.4.1)
##     Matrix               * 1.4-0      2021-12-08 [1] CRAN (R 4.2.2)
##     MatrixGenerics         1.18.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     matrixStats            1.4.1      2024-09-08 [1] RSPM (R 4.4.0)
##     memoise                2.0.1      2021-11-26 [1] RSPM
##     mgcv                   1.9-1      2023-12-21 [3] CRAN (R 4.4.1)
##     mime                   0.12       2021-09-28 [1] RSPM
##     miniUI                 0.1.1.1    2018-05-18 [1] RSPM (R 4.1.0)
##     modelr                 0.1.11     2023-03-22 [1] RSPM (R 4.4.0)
##     msigdbr              * 7.5.1      2022-03-30 [1] RSPM (R 4.1.0)
##     munsell                0.5.0      2018-06-12 [1] RSPM
##     nlme                   3.1-164    2023-11-27 [3] CRAN (R 4.4.1)
##     paletteer              1.6.0      2024-01-21 [1] RSPM (R 4.4.0)
##     parallelly             1.31.0     2022-04-07 [1] RSPM (R 4.1.0)
##     parameters             0.23.0     2024-10-18 [1] RSPM (R 4.4.0)
##     patchwork              1.3.0      2024-09-16 [1] RSPM (R 4.4.0)
##     pbapply                1.5-0      2021-09-16 [1] RSPM (R 4.1.0)
##     pillar                 1.9.0      2023-03-22 [1] RSPM (R 4.4.0)
##     pkgbuild               1.4.6      2025-01-16 [1] RSPM (R 4.4.0)
##     pkgconfig              2.0.3      2019-09-22 [1] RSPM
##     pkgload                1.3.4      2024-01-16 [1] RSPM (R 4.4.0)
##     plotly                 4.10.0     2021-10-09 [1] RSPM (R 4.1.0)
##     plyr                 * 1.8.7      2022-03-24 [1] RSPM (R 4.1.0)
##     png                    0.1-7      2013-12-03 [1] RSPM
##     polyclip               1.10-0     2019-03-14 [1] RSPM
##     prettyunits            1.1.1      2020-01-24 [1] RSPM (R 4.1.0)
##     print.rmd.tabs       * 0.2.0      2024-06-20 [1] Github (avishai987/print.rmd.tabs@ed6b345)
##     prismatic              1.1.0      2021-10-17 [1] RSPM (R 4.1.0)
##     processx               3.8.2.9000 2023-08-28 [1] Github (r-pkgs/processx@25812a1)
##     progress               1.2.3      2023-12-06 [1] RSPM (R 4.4.0)
##     promises               1.2.0.1    2021-02-11 [1] RSPM (R 4.1.0)
##     ps                     1.9.0      2025-02-18 [1] RSPM (R 4.4.0)
##     purrr                * 1.0.2      2023-08-10 [1] RSPM (R 4.4.0)
##     qvalue                 2.26.0     2021-10-26 [1] Bioconductor
##     R6                     2.5.1      2021-08-19 [1] RSPM
##     RANN                   2.6.1      2019-01-08 [1] RSPM (R 4.1.0)
##     rappdirs               0.3.3      2021-01-31 [1] RSPM
##     RColorBrewer           1.1-3      2022-04-03 [1] RSPM
##     Rcpp                   1.0.8.3    2022-03-17 [1] RSPM
##     RcppAnnoy              0.0.19     2021-07-30 [1] RSPM (R 4.1.0)
##     RCurl                * 1.98-1.6   2022-02-08 [1] RSPM
##     readr                * 2.1.5      2024-01-10 [1] RSPM (R 4.4.0)
##     readxl               * 1.4.3      2023-07-06 [1] RSPM (R 4.4.0)
##     rematch2               2.1.2      2020-05-01 [1] RSPM (R 4.1.0)
##     remotes                2.4.2      2021-11-30 [1] RSPM (R 4.1.0)
##     reprex                 2.1.1      2024-07-06 [1] RSPM (R 4.4.0)
##     reshape2               1.4.4      2020-04-09 [1] RSPM (R 4.1.0)
##     reticulate             1.42.0     2025-03-25 [1] RSPM (R 4.4.0)
##     rjson                  0.2.23     2024-09-16 [1] RSPM (R 4.4.0)
##     rlang                * 1.1.5      2025-01-17 [1] RSPM (R 4.4.0)
##     rmarkdown            * 2.29       2024-11-04 [1] RSPM (R 4.4.0)
##     ROCR                   1.0-11     2020-05-02 [1] RSPM (R 4.1.0)
##     roxygen2               7.1.2      2021-09-08 [1] RSPM (R 4.1.0)
##     rpart                  4.1.23     2023-12-05 [3] CRAN (R 4.4.1)
##     RSQLite                2.2.12     2022-04-02 [1] RSPM
##     rstatix              * 0.7.0      2021-02-13 [1] RSPM (R 4.1.0)
##     rstudioapi             0.17.1     2024-10-22 [1] RSPM (R 4.4.0)
##     Rtsne                  0.16       2022-04-17 [1] RSPM (R 4.1.0)
##     rvest                  1.0.4      2024-02-12 [1] RSPM (R 4.4.0)
##     S4Vectors            * 0.44.0     2024-10-29 [1] Bioconductor 3.20 (R 4.4.1)
##     sass                   0.4.10     2025-04-11 [1] RSPM (R 4.4.0)
##     scales                 1.3.0      2023-11-28 [1] RSPM (R 4.4.0)
##     scattermore            0.8        2022-02-14 [1] RSPM (R 4.1.0)
##     scatterpie             0.1.7      2021-08-20 [1] RSPM (R 4.1.0)
##     sctransform            0.3.3      2022-01-13 [1] RSPM (R 4.1.0)
##     sessioninfo            1.2.2      2021-12-06 [1] RSPM (R 4.1.0)
##     Seurat               * 4.1.0      2022-01-14 [1] RSPM (R 4.1.0)
##     SeuratObject         * 4.0.4      2021-11-23 [1] RSPM (R 4.1.0)
##     shadowtext             0.1.1      2022-01-10 [1] RSPM (R 4.1.0)
##     shape                  1.4.6.1    2024-02-23 [1] RSPM (R 4.4.0)
##     shiny                  1.9.1      2024-08-01 [1] RSPM (R 4.4.0)
##     SourceFromGithub     * 0.1.0      2022-12-07 [1] Github (avishai987/SourceFromGithub@99e0840)
##     spatstat.core          2.4-2      2022-04-01 [1] RSPM (R 4.1.0)
##     spatstat.data          3.0-1      2024-10-15 [1] Github (spatstat/spatstat.data@ae01346)
##     spatstat.geom          3.2-5      2024-10-15 [1] Github (spatstat/spatstat.geom@41263ad)
##     spatstat.random        3.1-5      2024-10-15 [1] Github (spatstat/spatstat.random@9528a39)
##     spatstat.sparse        3.0-1      2024-10-15 [1] Github (spatstat/spatstat.sparse@5fa3d52)
##     spatstat.utils         3.0-3      2024-10-15 [1] Github (spatstat/spatstat.utils@4e97cc3)
##     statsExpressions       1.6.1      2024-10-31 [1] RSPM (R 4.4.0)
##  VP stringi              * 1.7.12     2024-05-06 [?] RSPM (R 4.4.0) (on disk 1.8.4)
##     stringr              * 1.5.1      2023-11-14 [1] RSPM (R 4.4.0)
##     SummarizedExperiment   1.28.0     2022-11-01 [1] Bioconductor
##     survival               3.6-4      2024-04-24 [3] CRAN (R 4.4.1)
##     survminer              0.4.9      2021-03-09 [1] RSPM (R 4.1.0)
##     survMisc               0.5.6      2022-04-07 [1] RSPM (R 4.1.0)
##     svglite                2.1.3      2023-12-08 [1] RSPM (R 4.4.0)
##     systemfonts            1.1.0      2024-05-15 [1] RSPM (R 4.4.0)
##     TCGAbiolinks         * 2.35.3     2025-03-04 [1] Github (BioinformaticsFMRP/TCGAbiolinks@4075235)
##     TCGAbiolinksGUI.data   1.26.0     2024-10-31 [1] Bioconductor 3.20 (R 4.4.1)
##     tensor                 1.5        2012-05-05 [1] RSPM (R 4.1.0)
##     tibble               * 3.2.1      2023-03-20 [1] RSPM (R 4.4.0)
##     tidygraph              1.2.1      2022-04-05 [1] RSPM (R 4.1.0)
##     tidyr                * 1.3.1      2024-01-24 [1] RSPM (R 4.4.0)
##     tidyselect             1.2.1      2024-03-11 [1] RSPM (R 4.4.0)
##     tidytree               0.3.9      2022-03-04 [1] RSPM (R 4.1.0)
##     tidyverse            * 1.3.2      2022-07-18 [1] RSPM (R 4.4.1)
##     timechange             0.3.0      2024-01-18 [1] RSPM (R 4.4.0)
##     treeio                 1.18.1     2021-11-14 [1] Bioconductor
##     tweenr                 1.0.2      2021-03-23 [1] RSPM
##     tzdb                   0.5.0      2025-03-15 [1] RSPM (R 4.4.0)
##     ulimit                 0.0-3      2022-12-11 [1] Github (krlmlr/ulimit@67199f2)
##  VP units                * 0.8-1      2023-11-28 [?] RSPM (R 4.4.0) (on disk 0.8.5)
##     usethis              * 2.1.5      2021-12-09 [1] RSPM (R 4.1.0)
##     utf8                   1.2.2      2021-07-24 [1] RSPM
##     uwot                   0.1.11     2021-12-02 [1] RSPM (R 4.1.0)
##     vctrs                  0.6.5      2023-12-01 [1] RSPM (R 4.4.0)
##     viridis                0.6.2      2021-10-13 [1] RSPM
##     viridisLite            0.4.0      2021-04-13 [1] RSPM
##     vroom                  1.6.5      2023-12-05 [1] RSPM (R 4.4.0)
##     webshot                0.5.3      2022-04-14 [1] RSPM (R 4.1.0)
##     websocket              1.4.2      2024-07-22 [1] RSPM (R 4.4.0)
##     withr                  3.0.2      2024-10-28 [1] RSPM (R 4.4.0)
##     xfun                   0.52       2025-04-02 [1] RSPM (R 4.4.0)
##     XML                  * 3.99-0.9   2022-02-24 [1] RSPM
##     xml2                   1.3.3      2021-11-30 [1] RSPM (R 4.1.0)
##     xtable                 1.8-4      2019-04-21 [1] RSPM (R 4.1.0)
##     XVector                0.34.0     2021-10-26 [1] Bioconductor
##     yaml                   2.3.5      2022-02-21 [1] RSPM
##     yulab.utils            0.0.4      2021-10-09 [1] RSPM (R 4.1.0)
##     zeallot                0.1.0      2018-01-28 [1] RSPM (R 4.1.0)
##     zlibbioc               1.40.0     2021-10-26 [1] Bioconductor
##     zoo                    1.8-10     2022-04-15 [1] RSPM (R 4.1.0)
## 
##  [1] /sci/labs/yotamd/lab_share/avishai.wizel/Projects/libs
##  [2] /usr/local/lib/R/site-library
##  [3] /usr/local/lib/R/library
##  [4] /sci/home/avishaiw/R/x86_64-pc-linux-gnu-library/4.4
## 
##  V ── Loaded and on-disk version mismatch.
##  P ── Loaded and on-disk path mismatch.
## 
## ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

OPSCC_TCGA

Avishai Wizel

2025-09-02 15:07:00.929349

1 Data

2 MYB

3 HPV

4 MYB - HPV PPM correlation

5 HPV signature - HPV correlation

6 HPV signature - HPV T test

7 MYB- HPV T test

8 MYB and HPV survival plot

8.1 HPV + tumors

8.2 HPV- tumors

8.3 combine MYB and HPV survival plot

9 Myb Targets

10 HMSC signature

all samples

HPV+

HPV-

10.1 No cell cycle genes

OPSCC samples

HPV+

HPV-

10.2 T test

11 E2F ~ MYB

12 MYB ~ RBL1

13 OPSCC sc signature

13.1 all samples

13.2 HPV +

13.3 HPV -

14 Additional mutations

14.1 TP53

14.2 Fisher

14.3 TP53 unmutated

14.4 TP53 mutated

15 session info